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Abstract: We study the QCD phase diagram by first principle lattice calculations at so 
far unreached high densities. For this purpose we employ the density of states method. 
Unimproved staggered fermions, which describe four quark flavors in the continuum are 
used in this analysis. Though the method is quite expensive, small lattices show an indi- 
cation for a triple-point connecting three different phases on the phase diagram. 
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1. Introduction 

To clarify the phase diagram of QCD and thus the nature of matter under extreme condi- 
tions is one of the most interesting and fundamental tasks of high energy physics. The ex- 
perimental setup of the running Relativistic Heavy Ion Collider (RHIC) at the Brookhaven 
National Laboratory, as well as the forthcoming Large Hadron Collider (LHC) at CERN, 
aim on the creation of hot quark matter, the quark gluon plasma. Future experiments like 
the Compressed Baryonic Matter (CBM) experiment at the Facility for Antiproton and Ion 
Research (FAIR), attempt to explore the phase diagram in the region of high temperatures 
and intermediate densities. On the theoretical side various color-superconducting phases 
have been proposed in the region of very high densities but low temperatures, which may 
be relevant to the physics of neutron star interiors. For a review see for example [1]. 

Lattice QCD has been shown to provide important and reliable information from first 
principles on QCD at zero density. However, lattice QCD at finite densities has been 
harmed by the complex action problem ever since its inception. At non-zero chemical 
potentials (/x > 0) the determinant of the fermion matrix (detM) becomes complex. Stan- 
dard Monte Carlo techniques using importance sampling are thus no longer applicable 
when calculating observables in the grand canonical ensemble according to the partition 
function 

Z GC {fi)= / WdetM[U](/i)exp{-S G [U]}. (1.1) 

Note that Zqc is real. The origin of the sign problem are the fluctuations of the complex 
phase 9, defined by detM = |detM| exp{i#}. In case those fluctuations ( y ((9 — (9)) 2 )) 
become considerably larger than it/ 2 [2], the problem becomes serious. For a detailed 



discussion of the phase of the fermion determinant see also [3]. Recently many different 
methods have been developed to circumvent the complex action problem for p/T < 1 [4, 5]. 
For a recent overview see also [6]. With all these methods the transition line of the QCD 
phase transition was mapped in the (T,/z)-plane. The results coincide for small chemical 
potentials and the same choice of parameters. 

The goal of this paper is to employ the density of states method at non-zero chemical 
potential to extend the accessible region of the QCD phase diagram to lower temperatures 
and higher densities. The paper is organized as follows: In Section 2 we introduce our 
method and in Section 3 we give the details of our simulation parameters. Readers who 
are not interested in the details of our method or the simulation parameters may start at 
Section 4, where we will discussed results on the phase diagram. Results for the quark 
number density will be given in Section 5. Finally we will conclude in Section 6. 



2. Formulation of the method 

A very general formulation of the DOS method is the following: One special parameter ((f)) 
is held fixed. The expectation value of a thermodynamic observable (O), according to the 
usual grand canonical partition function (1.1), can be recovered by the integral over 0, 

<0>= J dcp (OfiU)}^^) I J dcf> </(£/))^) (2.1) 

where the density of states (p) is given by the constrained partition function: 

p(x) = Z^x) = j VUg(U) S(4> - x). (2.2) 

With ( we denote the expectation value with respect to the constrained partition func- 
tion. In addition, the product of the weight functions /, g has to give the correct measure 
of Zqc'- fg = detMexp{— Sq}- This idea of reordering the partition functions was used in 
many different cases [7, 8, 10, 11]. The advantages of this additional integration becomes 
clear, when choosing <j> = P and g(U) = 1. Here P denotes the plaquette expectation value. 
In this case p((f>) is independent of all simulation parameters. Once p(4>), as well as the 
correlation between <fi and the observable O is known, O can be calculated as a continuous 
function of the lattice coupling (3. If one has stored all eigenvalues of the fermion matrix 
for all configurations, the observable can also be calculated as a function of quark mass 
(m) and number of flavors [8] (Nf). 
In this work we chose 

(/) = P and g = \detM\ exp{-S G }, / = exp{i6}. (2.3) 

In other words we constrain the plaquette and perform simulations with measure g. This 
particular choice for the functions /, g was first introduced in ref. [9] and is called factoriza- 
tion method. It has been successfully tested on ^-vacuum like systems [10] and in random 
matrix theory [11]. By this choice, one has a clear separation of the effects coming from 
the complex phase of the determinant, and the rest. As we will see later, this enables us 



to locate the region of the parameter space which show the fewest fluctuation of the phase 
and hence contribute most to the partition function. 

In practice, we replace the delta function in eq. (2.2) by a sharply peaked potential 
[11]. The constrained partition function for fixed values of the plaquette expectation value 
can then be written as 

p(x) « JVU g(U) exp {-V(x)} , (2.4) 
where V(x) is a Gaussian potential with 

V{x) = ^{x-P) 2 . (2.5) 

We obtain the density of states (p(x)) by the fluctuations of the actual plaquette P around 
the constraint value x. The fluctuation dissipation theorem gives 

^lnp(x) =< 7 (x-P) > x . (2.6) 

Before performing the integrals in eq. (2.1) we carry out a calculation based on an ensemble 
generated at (ho,Po) : 

(Of(U)) x (h, f3) = (Of(U)R(p, no, f3, f3 )) x / (R(h, ho, P, Po)) x , (2.7) 
(f(U)) x (h, (3) = (f(U)R(ti, no, (3, (3o)) x I (R(h, ho, P, Po)) x , (2.8) 

^ In p(x, (3) = ( 7 (x - P)R{h, ho, P, Po)) x . (2.9) 

Here R is given by the ratio of the measure g at the point (h,/3) and at the simulation 
point (ho,Po), 

R(H,Ho,P,Po)=9(»,P)/g(»o,Po) = jp^exp{-S G ([3) + S G (Po)}- (2.10) 

|det(// )| 

Having calculated the expressions (2.7)-(2.9), we are able to extrapolate the expectation 
value of the observable (2.1) to any point (fi,(3) in a small region around the simulation 
point (ho,Po)- F° r an y evaluation of (0)(/z,/3), we numerically perform the integrals 
in eq. (2.1). We also combine the data from several simulation points as described by 
Eqs. (2.7)-(2.9). 

2.1 Simulations with constrained Plaquette 

The quantity we want to constrain is the real part of the plaquette P^uiy) averaged over 
lattice points (y) and directions (/x, v) 

P = E E ^TrP^(y). (2.11) 

V l<fJ.<u<4 



Since the plaquette is also the main part of the gauge action, 



S G = -(3J2 E {^ReTrP^(x)-l|, 



(2.12) 



x 1<^<u<4 



the additional potential V, which constraints the plaquette around a given value, can be 
easily introduced in the hybrid Monte Carlo update procedure of the hybrid-R algorithm 
[12]. To do so we modify the force in the molecular dynamical evolution of the gauge field 
by a factor (1 + j(x — P)/(3). This requires the measurement of the plaquette in each 
molecular dynamical step. 

2.2 Generating Configurations with measure |detM( / u)| 

For the generation of Configurations with measure g, according to eq. (2.3) we use the 
method of ref. [13]. Since at finite iso-spin chemical potential the fermion determinant is 
real, our fermion matrix in flavor space is given by 



Here each component represents a usual staggered fermion field \ with four flavors in the 
continuum limit. The diagonal elements are thus the usual staggered fermion matrices, in 
the upper left corner with chemical potential [i and in the lower right corner with chemical 
potential —fi. The off-diagonal elements are iso-spin symmetry breaking terms, propor- 
tional to the small parameter A, which are necessary in order to see spontaneous symmetry 
braking on a finite lattice. The fermion matrix (2.13) will thus represent eight continuum 
flavor. Note that the 75 matrices corresponds in the staggered case to a multiplication of 
the phase e(x) = -pi+z2+x 3 +s 4 . 

In order to simulate this system, we use the HMC in complete analogy to the even-odd 
ordering of the \i = case. This means that we generate Gaussian noise (R) for both 
components of {R^R^) = M(xi,X2), but keep only the upper component (xi) after the 
inversion. This way we still describe eight flavors with the block-diagonal positive definite 
matrix M^M which we use for the evolution of the gauge field. In the molecular dynamical 
integration of the equations of motion, we perform the usual square-root trick to reduce the 
number of flavors to four. In the limit of A — > it is however correct to take the square-root 
and will not introduce any approximations or locality problems, since our basic building 
block is the four flavor staggered fermion matrix which comes out with a power of one. A 
further reduction of rif to two or one, by using a fractional power of the staggered fermion 
matrix, introduces additional difficulties that are still under debate [14]. It is known [15] 
that using the square (or forth) root of the staggered matrix at \i > could lead to phase 
ambiguities . 

3. Simulation details and the strength of the sign problem 

Simulations have been performed with staggered fermions and Nf = 4. We chose 9 different 
points in the ((3, a / u)-plane for the 4 4 lattice and 8 points for the 6 4 lattice. On each of 
these points we did simulations with 20-40 constrained plaquette values, all with quark 
mass am = 0.05. Further simulations have been done with (f3,a[i) = (5.1,0.3) on the 
6 3 x 8 lattice for am = 0.05 and am = 0.03. The simulation points and statistics are 
summarized in table 1. For each point (/?, a/i) we have chosen several plaquette values in 
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Table 1: Simulation parameter and statistics. The number of trajectories is given in the last 
column, here the first factor is due to the number of the lambda parameters and the second factor 
gives the number of plaquette values. 

the range where the density of states (p) is large. In addition we measure all eigenvalues of 
the reduced fermion matrix. Here "reduced" means that the p dependence was reduced by 
Gauss elimination to the first and last time slice. This procedure was described in detail in 
ref. [4]. Having stored the eigenvalues, we are able to evaluate the fermionic determinant 
(absolute value and phase) as a function p, as well as all their derivatives. 

In order to calculate the plaquette expectation value, or its susceptibility, one has to 
perform the following integrals: 

(P) = j dx xp(x) (cos(d)) x , (P 2 ) = J dx x 2 p(x) (cos(0)) x . (3.1) 

Thus the functions p(x) and (cos(0)) x have to be known quite precisely. We plot both 
functions in figure 1. The transition is signaled in the double peak structure of p(x). 
The phase factor (cos(0)) x suppresses the peak of p(x) at smaller plaquette values, which 
results in a shift of the critical temperature to smaller values, in comparison with the phase 
quenched theory. 

In figure 2 we show the phase factor for different chemical potentials. With increasing 
chemical potential the phase factor becomes compatible with zero within errors. In fact, 
its average value becomes as low as cos(#) ~ 0.005. There exists however a small interval 
around P ~ 2.85, where the phase factor remains finite. In this way, the Plaquette ex- 




Figure 1: Results for simulations at f3 = 4.98, /i = 0.3, A = 0.02, rif = 4, am = 0.05, and number 
of lattice points: 4 4 . Shown are the density of states p(x), the phase factor (cos(0)), and their 
product. 
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Figure 2: Results for simulations at /? = 5.1, A = 0.01, n/ = 4, arri = 0.05, and number of lattice 
points: 6 4 . Shown is the suppression from the complex phase of the fermion determinant (cos(6)) 
for different chemical potentials. 



pectation values are strongly altered by the phase factor, figure 2 demonstrates also the 
advantage of he DOS method over the other approaches of lattice QCD at finite density. 
Using the DOS method one is able to do simulations directly at those Plaquette values 
which are relevant at finite density. This results in an even better overlap than that of the 
multi-parameter re-weighting approach. 
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Table 2: Results for the Sommer radius (ro), string tension (a), pion mass (m w ) and nucleon mass 
(mjv) for different (3 values. The results are for rif = 4 and am = 0.05, measured on a 10 3 x 20 
lattice. 



In order to be able to transform the lattice units into physical units, we have measured 
the heavy quark potential on a zero temperature lattice: 10 3 x 20. From the potential we 
determined the Sommer scale ro and the string tension a. In Addition we have measured 
the pion mass and the nucleon mass. The results are given in table 2. To transform 
the results from lattice to physical units in practice, the Sommer parameter ro was used 
and has been fitted with a polynomial of the order three, to get the lattice spacing a as 
a continuous function of (3. For ro in physical units we have used the MILC result of 
r = 0.476(7)(18) fm [16]. 

First of all we check, whether we can reproduce old results with our new method. For 
this purpose we reproduce one point of the phase transition line which has been calculated 
by multi-parameter re- weighting from [i = configurations on the 4 4 lattice [4] . Performing 
the integrations in eq. (3.1) numerically, we calculate the plaquette expectation value and 
the plaquette susceptibility xp = {P 2 ) ~ (P) 2 ■ At an y fixed A, we determine the critical 
coupling by the peak position of the susceptibility as shown in figure 3. We indeed find 
that including the phase factor does shift the transition to lower values of the coupling, 
which also means to lower temperatures. This can be clearly seen in a shift of the peak of 
the plaquette susceptibility. 

Since the A parameter introduces a systematic error, we have to perform a linear 
A — > extrapolation as shown in figure 4. The extrapolated result [3 = 4.938(4) (including 
the phase factor) and the result from multi-parameter re-weighting [4] are in very good 
agreement. The A dependence is expected to be smaller for larger /i, therefore from now 
on we only give results for X/m = 0.2. 
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Figure 3: Simulation parameters as in figure 1. Shown is the plaquette susceptibility as function 
of the coupling (3. 
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Figure 4: The A — > extrapolation of the critical couplings on the 4 4 lattice at a/i = 0.3. 



4. The Plaquette expectation value and the phase diagram 

Let us now discuss results for the plaquette expectation values from the 6 lattice as shown 
in figure 5 and figure 6. 

At chemical potentials \x < 0.36, the plaquette signals the QCD transition through 
a rapid crossover from a low temperature phase of < P >~ 2.9 to a high temperature 
phase of < P >~ 3.1. For \x > 0.36 the plaquette expectation value at small temperatures 
drops to < P >~ 2.85. This new low temperature phase of the plaquette at high chemical 
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Figure 5: Results for simulations at A = 0.01, nf =4, am = 0.05, and number of lattice points: 
6 4 . Shown is the Plaquette expectation value as a function of the coupling j3 for different chemical 
potentials 
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Figure 6: Simulation parameters as in figure 5. Shown is the plaquette expectation value at fixed 
coupling, as a function of the chemical potential. 



potentials is caused by the fermion determinant. As one can see in figure 2 the region 
around P ~ 2.85 is the region which is less suppressed by the phase factor. Another 
interesting observation is that the critical coupling, which is decreasing in \i for \i < 0.36 
starts to increase for \i > 0.36. 

The plaquette expectation value thus suggests the existence of three different phases in 
the (T,/i)-diagram with a triple point, where all those phases coincide. In figure 7 we show 
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Figure 7: The phase diagram in physical units from Nt = 6 and 8 lattices. 

the phase diagram in physical units. The phase boundaries were determined by calculating 
the peaks in the plaquette susceptibility. The points at T = 93 MeV are from calculations 
on 6 3 x 8 lattices. Note, that we make no statement about the order of the transition lines. 
To determine the order of the transition one has to perform a finite-size-scaling analysis 
which is beyond the scope of this article. 

The triple point is located around /j,^ 1 300 MeV, however its temperature (T tri ) 
decreases from T tri 148 MeV on the 4 4 lattice to T tri « 137 MeV on the 6 4 lattice. This 
shift reflects the relatively large cut-off effects one faces, with standard staggered fermions 
and temporal extents of 4 and 6. 

Also shown in figure 7 are points from simulations with quark mass am = 0.03. The 
phase boundary at low temperatures turned out to be — within our statistical uncertainties 
- independent of the mass. This gives evidence that the transition is associated with 
the onset of baryonic matter rather than a pion condensate. Going from am = 0.05 to 
am = 0.03, the pion mass changes by a factor of ~ y/3/5 = 0.77 whereas the nucleon mass 
remains approximately constant. 
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5. The quark number density 

To reveal the properties of the new phase located in the lower right corner of the phase 
diagram, we calculated the quark number density, at constant coupling (3 and at constant 
temperature respectively. To obtain the density n q we perform the following integration 



dx (— - — — cos(0)> p(x) (5.1) 
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Figure 8: The quark number density at constant temperature T = 143 MeV (4 4 lattice), T 
124 MeV (6 4 lattice) and T = 93 MeV (6 3 x 8 lattice). 



The thermodynamic quantity n q are given as usual by 

1 /dlndetM 
Uq ~ a 3 NiN t \ d(a^) 

In figure 8 we show the baryon number density, which is related to the quark number density 
by ns = n q /3. The results are plotted in physical units and correspond to a constant 
temperature of T « 143 MeV (4 4 lattice), T w 124 MeV (6 4 lattice) and T sa 93 MeV 
(6 x 8 lattice). In order to remove the leading order cut-off effect, we have multiplied the 
data with the factor c = SB(Nt)/SB, which is the Stefan-Boltzmann value of a free lattice 
gas of quarks at a given value of Nt, divided by its continuum Stefan-Boltzmann value. 
For unimproved staggered fermions this correction factor can be as large as 2 and will not 
completely remove all the cutoff effects. 

At the same value of the chemical potential where we find also a peak in the suscepti- 
bility of the plaquette (/U c ), we see a sudden rise in the baryon number density. Thus for 
fi>fj> c we enter a phase of dense matter. The transition occurs at a density of (2 — 3) x njv, 
where denotes nuclear matter density. Above the transition, the density reaches values 
of (10 — 20) x Tijy. Quite similar results have been obtained recently by simulations in the 
canonical ensemble [17]. 

6. Conclusions 

We have explored the QCD phase diagram by first principle lattice calculations at so far 
unreached hight densities. In the accessible region of T > 100 MeV and fx q < 400 MeV we 
have been able to identify three different regions, which seem to be separated by different 
plaquette expectation values and quark number densities. These regions coincide in a 




triple-point. The triple-point has been located at T tri < 137 MeV and m q ~ 300 MeV at 
finite lattice spacing (N t = 6). We note, that the lowest reachable temperature on N t = 6, 
is about 60 MeV (setting the scale e.g. by m p ). We thus find the triple-point about 80 MeV 
above the lowest temperature. This is the first numerical evidence from lattice QCD for 
a third phase (appart form the hadronic phase and the quark gluon plasma) and a triple 
point in the QCD phase diagram. In the third phase the quark number density reaches 
values of (10 — 20) x tin, where tin denotes normal nuclear matter density. 

The new phase is a natural candidate for a color superconducting phase. Recently, 
by combining experimental results from cold atoms in a trap [18] and some universal 
arguments, an upper bound for the transition line from the quark gluon plasma phase 
(QGP) to the superconducting phase (SC) was proposed (T c < 0.35Ep) [19, 20]. To first 
approximation the Fermi-Energy Ep is given by the chemical potential fi q . In [20] the 
triple point was estimated by comparing this upper bound with the experimental freeze- 
out curve. A temperature of T tn < 70 MeV was found. Our value of the triple-point 
roughtly corresponds to T c < OAGEp- It is interesting that the two values are close. 

At low temperatures we find a phase boundary which is very steep and almost inde- 
pendent of T. Although our lowest temperature is 96 MeV an extrapolation to T = seems 
to be reasonable and would yield a critical chemical potential of n q (T = 0) ~ 250 MeV or 
equivalently hb/T c (/j,b = 0) ~ 4.7. This number appears to be at the lower edge of the 
phenomenological expectation of hb/T c ([Xb = 0) ~ 5 — 10. Note, that our lattice spacing 
is close to the strong coupling regime and we should feel the influence of the strong cou- 
pling limit. Strong coupling expansion calculations in general yield much lower values of 
» b /T c <j*b = 0) < 1.5 [21]. 

For this work the density of state method has been employed, which works well on 
small lattices up to chemical potentials of n q /T < 3 (other methods [4, 5] worked up to 
Hq/T < 1). The method is however extremely expensive and thus will in the near future not 
yield results close to the thermodynamic limit or the continuum limit, due to limitations 
in computer resources. 

We have to emphasize once more, that the simulations have been carried out on coarse 
lattices with an unphysical value of rif = 4 degenerate fermion flavor, and that neither 
the continuum nor the thermodynamic limit has been taken. Since we used unimproved 
staggered fermions, the corrections due to a finite lattice spacing are large. We also expect 
corrections due to the finite size of our volume. The simulations have not been performed 
with a constant quark mass, but m q /T = 0.3 has been held fixed. 
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